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Longitudinal imaging studies are essential to understanding the 
neural development of neuropsychiatric disorders, substance use dis- 
' orders, and the normal brain. The main objective of this paper is 

to develop a two-stage adjusted exponentially tilted empirical like- 
ly ' lihood (TETEL) for the spatial analysis of neuroimaging data from 
, longitudinal studies. The TETEL method as a frequentist approach 

allows us to efhciently analyze longitudinal data without modeling 
temporal correlation and to classify different time-dependent covari- 
^ ■ ate types. To account for spatial dependence, the TETEL method 

' developed here specifically combines all the data in the closest neigh- 

, borhood of each voxel (or pixel) on a 3-dimensional (3D) volume (or 

CO ' 2D surface) with appropriate weights to calculate adaptive parame- 

(N ■ ter estimates and adaptive test statistics. Simulation studies are used 

QQ ■ to examine the finite sample performance of the adjusted exponen- 

c 2 ^ ' tial tilted likelihood ratio statistic and TETEL. We demonstrate the 

application of our statistical methods to the detection of the differ- 
ence in the morphological changes of the hippocampus across time 
between schizophrenia patients and healthy subjects in a longitudinal 
L J ■ schizophrenia study. 
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1. Introduction. Neuroimaging data, including both anatomical and 
functional magnetic resonance imaging (MRI), have been/are being widely 
collected to understand the neural development of neuropsychiatric disor- 
ders, substance use disorders, and the normal brain in various longitudinal 
studies [Almli et al. (2007)]. For instance, various morphometrical measures 
of the morphology of the cortical and subcortical structures (e.g., hippocam- 
pus) are extracted from anatomical MRIs for understanding neuroanatom- 
ical differences in brain structure across different populations and across 
time. Studies of brain morphology have been conducted widely to char- 
acterize differences in brain structure across groups of healthy individuals 
and persons with various diseases, and across time [Thompson and Toga 

(2002) , Thompson, Cannon and Toga (2002), Styner et al. (2005), Zhu et 
al. (2008a)]. Moreover, functional MRI (fMRI) is a valuable tool for un- 
derstanding functional integration of different brain regions in response to 
specific stimuli and behavioral tasks and detecting the association between 
brain function and covariates of interest, such as diagnosis, behavioral task, 
severity of disease, age, or IQ [Priston (2007), Rogers et al. (2007), Huettel, 
Song and McCarthy (2004)]. 

Much effort has been devoted to developing frequentist and Bayesian 
methods for analyzing neuroimaging data using numerical simulations and 
theoretical reasoning. Frequentist statistical methods for analyzing neu- 
roimaging data are often sequentially executed in two steps. The first step 
involves fitting a general linear model or a linear mixed model to neuroimag- 
ing data from all subjects at each voxel [Beckmann, Jenkinson and Smith 

(2003) , Friston et al. (2005), Rowe (2005), Woolrich et al. (2004), Zhu et 
al. (2008a)]. The second step is to calculate adjusted values that account 
for testing the hypotheses across multiple brain regions or across many vox- 
els of the imaging volume using various statistical methods (e.g., random 
field theory, false discovery rate, or permutation method) [Cao and Wors- 
ley (2001), Friston et al. (1996), Hayasaka et al. (2004), Logan and Rowe 

(2004) , Worsley et al. (2004)]. Most of these frequentist methods have been 
implemented in existing neuroimaging software platforms, including statis- 
tical parametric mapping (SPM) (www.fil.ion.ucl.ac.uk/spm/) and FMRIB 
Software Library (FSL) (www.fmrib.ox.ac.uk/fsl/), among many others. In 
the recent literature, a number of papers have been published on the de- 
velopment of Bayesian spatial-temporal models for functional imaging data 
[Penny, Flandin and Trujillo-Barreto (2007), Bowman et al. (2008), Woolrich 
et al. (2004), Luo and Puthusserypady (2005)]. Most Bayesian approaches, 
however, are less practical due to the extensively computational burden of 
running a Markov chain Monte Carlo method in a large number of vox- 
els [Bowman et al. (2008)], and, thus, they are limited to small or moderate 
anatomic regions and a small number of regions of interest (ROI). Moreover, 
as pointed out in Snook, Plewes and Beaulieu (2007), the major drawbacks 
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of ROI analysis include the instability of statistical results obtained from 
ROI analysis and the partial volume effect in relative large ROIs. 

Existing statistical methods in the neuroimaging literature have two ma- 
jor limitations for analyzing longitudinal neuroimaging data, as explained 
below. The respective strategies to resolve these two limitations are detailed 
in Section 2. The first limitation is that the parametric models, such as 
linear mixed models, require the correct specification of the temporal corre- 
lation structure and cannot properly distinguish between different types of 
time-dependent covariates (types I, II and III) [Diggle et al. (2002), Lai and 
Small (2007), Pepe and Anderson (1994)]. A distinctive feature of longitu- 
dinal neuroimaging data is that it is able to characterize individual change 
in neuroimaging measurements (e.g., volumetric and morphometric) over 
time, and the time-dependent covariates of interest may influence change. 
Imaging measurements of the same individual usually exhibit positive corre- 
lation and the strength of the correlation decreases with the time separation 
[Liang and Zeger (1986)]. Moreover, longitudinal data may provide crucial 
information for a causal role of a time-dependent covariate (e.g., exposure) 
in the disease process [Diggle et al. (2002), Lai and Small (2007), Pepe 
and Anderson (1994)]. Improperly handling time-dependent covariates and 
ignoring (or incorrectly modeling) temporal correlation structure in imag- 
ing measures likely would influence subsequent statistical inference, such as 
increasing the false positive and negative errors, and result in misleading 
scientific inferences [Diggle et al. (2002), Lai and Small (2007)]. 

The second limitation is that most smoothing methods apply the same 
amount of smoothing throughout the whole image, which can be problem- 
atic near the edges of the significant regions. Although it is common to 
apply a smoothing step before applying a voxel-wise approach for the anal- 
ysis of neuroimaging data [Poline and Mazoyer (1994), Shafie et al. (2003), 
Lindquist and Wager (2008)], the voxel- wise method suffers from the same 
amount of smoothing throughout the whole image and the arbitrary choice 
of smoothing extent [Hecke et al. (2009), Jones et al. (2005)]. Jones et al. 
(2005) have shown that the final results of voxel-based analysis can strongly 
depend on the amount of smoothing in the smoothed diffusion imaging data. 
Recently, Yue, Loh and Lindquist (2010) introduced a spatially smooth- 
ing method using nonstationary spatial Gaussian Markov random fields to 
spatially and adaptively smooth images. Their approach, however, can be 
computationally extensive for 3D imaging data. 

In this paper we will develop strategies to resolve these two limitations. 
To resolve the first limitation, we develop an adjusted exponentially tilted 
empirical likelihood method, called AETEL, for the analysis of longitudinal 
neuroimaging data with time-dependent covariates. AETEL is a nonpara- 
metric method that is built on a set of estimating equations and the number 
of estimating equations can be larger than the number of parameters. Thus, 
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it avoids parametric assumptions and this feature is very appealing for the 
analysis of real neuroimaging data, such as brain morphological measures, 
because the distribution of the univariate (or multivariate) neuroimaging 
measurements often deviates from the Gaussian distribution [Ashburner and 
Friston (2000), Salmond et al. (2002), Luo and Nichols (2003)]. Using more 
estimating equations than the number of parameters allows us to appropri- 
ately handle time-dependent covariates of different types and to make an 
efficient use of the estimating equations without the need of modeling the 
temporal correlation in longitudinal data [Lai and Small (2007), Qu, Lindsay 
and Li (2000)]. AETEL also provides a natural test statistic to test whether 
a specific covariate is of a certain type (types I, II and III). 

To resolve the second limitation, we develop a two-stage AETEL, abbrevi- 
ated as TETEL, for the analysis of longitudinal neuroimaging data. TETEL 
integrates a smoothing method into our AETEL for carrying out statisti- 
cal inference on neuroimaging data. The TETEL method, as an adaptive 
procedure, fits AETEL at each voxel in stage 1. Then, TETEL uses the 
information learned from stage 1 to discard the data from the neighboring 
voxels with dissimilar signal pattern and to incorporate the data from the 
neighboring voxels with similar signal pattern to adaptively calculate param- 
eter estimates and test statistics. TETEL avoids using the same amount of 
smoothing throughout the whole image in most smoothing methods. In ad- 
dition, theoretically, we can establish consistency and asymptotic normality 
of the estimators and test statistics obtained from TETEL. 

Section 2 of this paper introduces the shape data of the hippocampus 
structure from a longitudinal schizophrenia study and presents the new sta- 
tistical methods just described. In Section 3 we conduct simulation stud- 
ies to examine the finite sample performance of the TETEL method. Sec- 
tion 4 illustrates an application of the proposed methods to the longitudinal 
schizophrenia study of the hippocampus. We present concluding remarks in 
Section 5. 

2. Data and methods. 

2.1. Longitudinal schizophrenia study of hippocampus shape. This is a lon- 
gitudinal, randomized, controlled, multisite, double-blind study conducted 
at 14 academic medical centers in North America and western Europe, with 
partial funding from Lilly Research Laboratories [Lieberman et al. (2005), 
Styner et al. (2004)]. In this study 238 first-episode schizophrenia patients 
were enrolled meeting the following criteria: age 16 to 40 years; onset of psy- 
chiatric symptoms before age 35; diagnosis of schizophrenia, schizophreni- 
form, or schizoaffective disorder according to the fourth edition of diagnostic 
and statistical manual of mental disorders (DSM-IV) criteria; and various 
treatment and substance dependence conditions. After random allocation at 
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Fig. 1. Location of hippocampus in the context of the surrounding structures in the coro- 
nal (a) and sagittal (b) views. Subregions of the hippocampus in (c) showing the head of 
the hippocampus (HH), the digitationes hippocampi (DH), the hippocampal body (HB), the 
hippocampal tail (HT), the terminal segment of the HT (TS), the dentate gyrus (DG), 
and the fields of the comu ammonis (CA1-CA4)- Adapted with permission from Springer 
Verlag, Heidelberg, Germany [Duvernoy (2005)]. 

baseline, 123 patients were selected to receive a conventional antipsychotic, 
haloperidol (2-20 mg/d), and 115 were selected to receive an atypical an- 
tipsychotic, olanzapine (5-20 mg/d). Patients were treated and followed up 
to 47 months. Also, 56 healthy control subjects matched to the patient's 
demographic characteristics were enrolled. Neurocognitive and MRI assess- 
ments were performed at months (baseline), 3, 6, 13, 24, 36, and 47 ap- 
proximately, with different subjects having different visiting times, and some 
subjects dropped out during the course of the study. 

The hippocampus, a gray matter structure in the limbic system, is in- 
volved in processes of motivation and emotions and has a central role in the 
formation of memory. The hippocampus is a paired structure with mirror- 
image halves in the left and right brain hemispheres and located inside the 
medial temporal lobe (Figure 1). Many MRI studies have reported the re- 
duction of hippocampal volume demonstrated in schizophrenia subjects and 
at onset of the first episode of psychotic symptoms before effects associated 
with treatment and disease chronicity [Lieberman et al. (2005)]. 

The aim of this study is to use the boundary and medial shape of the 
hippocampus to examine whether hippocampal abnormalities are present in 
schizophrenia patients. Statistical shape modeling and analysis have emerged 
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as important tools for understanding cortical and subcortical structures from 
medical images [Dryden and Mardia (1998)]. We consider two approaches 
for shape representation including a spherical harmonic description sampled 
into a triangulated surfaces (SPHARM-PDM) and a medial shape descrip- 
tion [Pizer et al. (2003), Styner and Gerig (2003)]. The SPHARM-PDM can 
only represent objects of spherical topology, whereas the medial representa- 
tion provides information on a rich set of features, including local thickness. 
These shape features are not accessible by conventional volume-based mor- 
phometry and offer us a great opportunity to address the weaknesses of 
conventional volumetric methods. 

We consider two sets of responses of interest. The first set of responses 
was based on the SPHARM-PDM representation of hippocampal surfaces. 
We use the SPHARM-PDM [Styner et al. (2004)] shape representation to es- 
tablish surface correspondence and align the surface location vectors across 
all subjects. The sampled SPHARM-PDM is a smooth, accurate, fine-scale 
shape representation (Figure 3). The hippocampal surfaces of different sub- 
jects are thus represented by the same number of location vectors (with each 
location vector consisting of the spatial x,y, and z coordinates of the corre- 
sponding vertex on the SPHARM-PDM surface) and are used as the second 
set of responses. Covariates of interest are race (Caucasian, African Ameri- 
can, and others), age (in years), gender, group (the schizophrenia group and 
the healthy control group) and time (visiting times in months). 

The second set of responses was the hippocampus m-rep thickness at the 
24 medial atoms of the left and the right brain (Figure 4). The m-rep is 
a linked set of medial primitives named medial atoms, which are formed 
from two equal length vectors and are composed of a position, a radius, 
a frame implying the tangent plane to the medial manifold, and an object 
angle [Styner et al. (2004)]. The m-rep thickness is the radius of each medial 
atom. Covariates of interest were WBV, race (Caucasian, African American, 
and others), age (in years), gender, diagnostic status (patient or control), and 
visiting times (in weeks). This WBV measure includes gray and white mat- 
ter, ventricular cerebrospinal fluid, cisterns, fissures, and cortical sulci. The 
WBV is commonly used as a covariate in statistical analyses to control for 
scaling effects [Arndt et al. (1991)]. Particularly, WBV is a time-dependent 
covariate and may vary with the hippocampus thickness measurement. 

2.2. Estimating equations for longitudinal data. We consider a longitudi- 
nal study of imaging data with n subjects, where a qxl covariate Xjj (e.g., 
age, gender, height, and brain volume) is obtained for the ith subject at the 
jth time point tij for i = 1, . . . ,n and j = 1,. . . , mi. Thus, there are at least 
Y17=i nT'i = ^ images in the study. Based on each image, we observe or com- 
pute neuroimaging measures, denoted by Yj = {yij{d) : d G'D,j = 1, . . . , rui}, 
across all rui time points from the ith subject, where d represents a voxel (or 
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atom, or point) on D, a specific brain region. The imaging measure yij{d) at 
each voxel d can be either univariate or multivariate. For example, the m-rep 
thickness is a univariate measure, whereas the location vector of SPHARM- 
PDM is a three-dimensional MRI measure at each point [Styner and Gerig 
(2003), Chung, Dalton and Davidson (2007)]. For notational simplicity, we 
assume that the yij{d) are univariate measures. 

We temporarily drop voxel d from our notation. At a specific voxel d in 
the brain region, Zj = {(yjj,Xjj) : j = 1, . . . ,mi} is independent and satisfies 
a moment condition 

(2.1) E{g{z„9)} = for i = l,...,n, 

where 9 is a p x 1 vector, g{-,-) is an r x 1 vector of known functions with 
r >p, and E denotes the expectation with respect to the true distribution of 
all the Zj's. Equation (2.1) is often referred to as a set of unbiased estimating 
equations or moments model [Qin and Lawless (1994), Hansen (1982)]. The 
moments model (2.1) is more general than most parametric models including 
linear mixed model used for the analysis of neuroimaging data [Worsley et 
al. (2004), Qin and Lawless (1994), Hansen (1982), Schennach (2007), Owen 
(2001)]. 

For longitudinal data, although the measurements from different sub- 
jects are independent, those within the same subject may be highly cor- 
related. The generalized estimating equations (GEE) assume a working co- 
variance matrix for = (y^i, . . . ,yim,V given by Vi. Let E{yi) = = 
(/^ii(/3), • ■ • , A*imi(/3))'^ and Di{f3) = d^i{f3)/d(3. Under the assumption that 
E{Di{l3)'^Vf^[yi — Aii(/3)]} = 0, Liang and Zeger (1986) proposed to use an 
estimator, denoted by /3gee, which solves a set of GEEs as follows: 

n 

(2.2) G(/3) = Y^D,i(3fV-'[y, - /z,(/3)] = 0. 

i=l 

For longitudinal data with time-dependent covariates, whether E[g(zi, 
9)] = E{Di{l3)'^V^^[yi — equals zero or not depends on the type of 

time-dependent covariates and the structure of Vi [Lai and Small (2007)]. 
The time-dependent covariate Xjj is of type I if 

(2.3) E{di3His{/3)[yij - fj,ij{(3)]} = for all s, j = 1, . . . , mj, 

where 9^ = d/dfi. A sufficient condition for type I covariates is £'[yjj|xjj] = 
E[yij\yiii,. . . ,XjmJ- For type I covariates, we can set g{'Zi,9)=Di{j3YV~^[yi- 
//j(/3)] and show that E[g{zi,9)\ = 0. If Vi is the covariance matrix of y,, 
then the estimator /3gee is an efficient estimator. However, /3gee is ineffi- 
cient under a misspecified Vi. To increase the efficiency, we may choose 
several candidate working covariance matrices mI^\. . . , ikf^^**"^ and assume 
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V^~^ = X^^^]^ UkM- ^^ for some unknown constants ak [Qu, Lindsay and Li 
(2000)]. Then, following Qu, Lindsay and Li (2000), we consider a set of 
estimating equations given by 



(2.4) g{z,,e) 



/A(/3)^Mf)[y,-/.,(/3)]\ 

VA(/5)^Mf°'^[y-/^.(/5)]/ 



f or i = 1 , 



,n. 



In this case, the number of functions in g{zi,6) is sqq > q, when sq > 1. 
The time-dependent covariate Xjj is of type II if 

(2.5) E{df}fj,is{(3)[yij - fiijiP)]} = for all s > j, j = 1, . . . , m^. 
A sufficient condition for type II covariates is 

(2.6) p{xi^t+i,-- ■ ,Xim^\yit,Xit) =p{xi^t+i, ■ ■ . ,XimJxit). 

For type II covariates, we can set g{zi,d) = A(/5)"^[yj — /ii (/?)], in which an 
independent working covariance matrix is used. However, the estimator /3gee 
based on the independent working correlation matrix is inefficient, since we 
do not use the information contained in E{dj3fiis{l3)[yij — fJ-ijif^)]} = for 
all s > j. To increase the efficiency of the estimate, we choose a set of lower 
triangular rrii x mj matrices L^^^ , . • • , , and then we consider estimating 
equations given by 



(2.7) 



for i = 1, . . . , 



n. 



In this case, the number of functions in g{zi,9) is sqq > q, when sq > 1. 
Supposing that mi = • • • = m„, we can set sq = mi (mi + l)/2 and L^-''^ = 
egej for s > j and 6 = 1, . . . , sq, where is a g x 1 vector with the sth 
component 1 and otherwise. Thus, similar to Lai and Small (2007), we are 
able to pick di3His{fi)[yij - ^J-iji(3)] for all s > j. 
The time-dependent covariate Xjj is of type III if 

(2.8) E{^p^iis{p)[y^j - fiijiP)]} / for some s > j. 

For type III covariates, we need to choose as a diagonal matrix. For in- 
stance, if Vi = Ij, where Ij is an x rrii identity matrix, then g{zi,9) = 
DiH^yiYi ~ Furthermore, if we assume the specific form for the vari- 

ances of all Yij, then we may set Vi = diag(Cov(yi)). 

An overall strategy to analyze models with time-dependent covariates is 
first to assume that the time-dependent covariates are of type III. Then we 
test whether the time-dependent covariates are of type II, and if the test 
is not rejected, we can go on to test if they are of type I. Once the type 
of all the time-dependent covariates is decided, we use the corresponding 
estimating equations. See Section 4 for more details. 
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2.3. Adjusted exponentially tilted empirical likelihood. We consider a non- 
parametric method, called an exponentially tilted empirical likelihood, to 
carry out statistical inference about 9 based on a set of estimating equations 
{g{zi,9) -.i = l,...,n} [Schennach (2007)]. The exponentially tilted empiri- 
cal likelihood (ETEL) method is a combination of the exponentially tilted 
(ET) method and the empirical likelihood (EL) method. Both EL [Owen 
(2001), Qin and Lawless (1994)] and ET [Imbens, Spady and Johnson (1998)] 
methods combine the reliability of nonparametric methods with the effec- 
tiveness of the likelihood approach. The EL estimator exhibits desirable 
higher-order asymptotic properties, whereas the EL estimator may fail to 
be -^/n-convergent in the presence of model misspcification. In contrast, the 
ETEL estimator maintains -^/n-convergence under model misspecification 
[Schennach (2007)]. 

However, most empirical likelihood type methods including ETEL suffer 
from two pitfalls: relatively low precision of the chi-square approximation 
and nonexistence of solutions to the estimating equations [Chen, Variyath 
and Abraham (2008), Liu and Chen (2010)]. Chen, Variyath and Abraham 
(2008) introduce a novel adjustment to these empirical likelihood meth- 
ods and develop an iterative algorithm that converges very fast. Simulation 
studies have shown that the adjusted empirical likelihood methods perform 
as well as the linear regression model with Gaussian noise when data are 
symmetrically distributed, while the adjusted empirical likelihood methods 
are superior when data have skewed distribution [Zhu et al. (2009), Chen, 
Variyath and Abraham (2008), Liu and Chen (2010)]. 

Following Chen, Variyath and Abraham (2008), we consider an adjust- 
ment of ETEL, abbreviated as AETEL, by introducing an adjustment 



n 



(2.9) 




i=l 



where an = max(l,log(n)/2). Then, AETEL is defined as 



n+l 



(2.10) 




i=l 



where Pi{9) is the solution to 



n+l 



min (n + 1) 
Pi,...,p„+i 




i=l 



subject to 



n+l 



n 




and 




i=l 
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The maximum AETEL estimator, denoted by ^Aeteb minimizes a criterion 
given by 

(2.11) ^Aetei = argmin£Aetei(6')- 

e 

According to a duality theorem in convex analysis [Newey and Smith (2004)], 

. eMKOYgn+m) exp(t(g)^ff(z,,g)) 
Pn^r{e) = and = 

for i = 1, . . . , n, in which 

TM = j;exp(t(0)^5(zi,^)) + exp(i(e)V+i(^)), 
i=i 

i{9) = argmax|-^exp(-t^5(zi,6i)) -exp(-t^5-„+i(6l))|. 

We use the numerical algorithm proposed by Chen, Variyath and Abraham 
(2008) to compute ^Aeteh which combines the modified Newton-Raphson al- 
gorithm and the simplex method. Compared with that of computing ETEL, 
this numerical algorithm of Chen, Variyath and Abraham (2008) converges 
very fast and the solution to AETEL is guaranteed. 
We consider testing the linear hypotheses: 

(2.12) Ho:Re = ho vs. HiiRe^ho, 

where i? is a cq x p matrix of full row rank and bo is a cq x 1 specified vector. 
Most scientific questions in neuroimaging studies can be formulated into 
linear hypotheses, such as a comparison of brain regions across diagnostic 
groups and a detection of changes in brain regions across time. The AETEL 
ratio statistic for testing R9 = bo can be constructed as follows: 

(2.13) MActel = -2(n+l)| sup £Aete\{0)-Snp£Actc\{0)}. 

Thus, to compute Li^Aetob we also need to compute the maximum AETEL 
estimator, denoted by ^Actci.Oi subject to an additional constraint RQ = \><^. 

Under some conditions on g{zi,9), we have the following theorem, whose 
detailed proof can be found in a supplementary document [Shi et al. (2011)]. 

Theorem 2.1. // assumptions (Al)-(A4) in the supplementary docu- 
ment are true, then we have the following: 

(a) \/n{9 Actci — (^o) converges to i'q = N{0,T,) in distribution, where 6q 
denotes the true value of 6 and T, = {DV~^D'^)^^ , in which 

n n 

D= lim n~^y"dgg{z„e) and V= lim V c/(zi, 0)®^; 

1=1 1=1 
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(b) under the null hypothesis Hq, LR\ctcl converges to a x^(co) distribu- 
tion; 

(c) if E[g{zi,9)] = for all i andr>p, then LRq-f =—2{n+l) supg iActci{(^) 
is asymptotically x^{r—p). 

We have established consistency and asymptotic normality of ^Aetei and 
the asymptotic distribution of LR^ctcX- Theorem 2.1 also shows that AE- 
TEL has the same first-order asymptotic properties as ETEL [Schennach 
(2007)]. High-order precision of AETEL can be explored by following the 
arguments in Liu and Chen (2010). It will be shown that the chi-square ap- 
proximation of the AETEL likelihood ratio statistics is found precise, com- 
pared with the existing ETEL [Owen (2001), Liu and Chen (2010), Chen, 
Variyath and Abraham (2008)]. Providing a reliable p- value at each voxel 
is crucial for controlling the family-wise error rate and false discovery rate 
(FDR) across the entire brain region [Benjamini and Hochberg (1995), Wors- 
ley et al. (2004)]. 

2.4. Two-stage adaptive estimation procedure. We propose a two-stage 
adaptive estimation procedure for computing parameter estimates and like- 
lihood ratio statistics for the spatial and adaptive analysis of neuroimaging 
data in 3D volumes (or 2D surfaces). To distinguish data and parameter 
in different voxels, we introduce voxel d into our notation. For instance, 
Zi{d) and 6{d), respectively, denote the zth observation and the parameter 
at voxel d. 

Stage 1 is to calculate the AETEL estimator of the parameter 0{d), de- 
noted by 9xete\{d), based on a set of estimating equations {g{zi{d),6{d)) : i = 
1, . . . , n} at each voxel dGD. 

One chooses a set of estimating equations {g{zi{d),9{d)) -.i = l,...,n} 
according to a specific type of time-dependent covariate and then substitutes 
them into (2.10) to build i\etei{6{d);d). Subsequently, we solve OActciid) 
according to (2.11) by minimizing £Actc\i(^id)',d), and then we obtain a set 

of parameter estimates {^Aetel(rf) : d € V}. 

Stage 2 is to calculate the TETEL estimator of 9{d), denoted by ^Tctci('^)) 
by utilizing the information contained in {OActciid) '■ d G T>} . Then, we cal- 
culate a TETEL ratio statistic, denoted by X-RTctei('^); for testing HQ[d): 
Re{d) = hQ. 

Specifically, one combines all data in the voxel d and the set of the closest 
neighboring voxels of d, denoted by N{d), to form a new set of estimating 
equations {g{zi{d),6{d);d) : i = 1, . . . , n} as follows: 

(2.14) ~giziid),e{d)-d)= oj{d';d)g{z,id'),e{d)), 

d'(^N{d)U{d} 
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where oj{d';d) is a weight describing the similarity between voxel d and 
any d' G N{d). The weights uj{d']d) at each d depend on the parameters 
{Oxete\{d') '■ d' G N{d) U {d}} calculated in Stage 1. From now on, we assume 
that u}{d';d) takes the form 

(2.15) uj{d'- d) = exp{-LRAetciid'; d)/Cn), 

where C„ = Xi_q,(p) log(n)/5 and Xi-aip) is the upper a-percentile of the 
X^ip) distribution. In addition, 

(2.16) LRAetci{d';d) = -2{n + l)\iAetei{OAetei{d');d)-snj>£Aetci{0;d)}, 

in which ^Actci(^if^) is defined in (2.10) based on the estimating equations 
{g{zi{d),6{d)) -.1 = 1,. . . ,n}. Statistically, LRj^ctciid'',d) denotes the AETEL 
ratio statistic for testing the hypothesis HQ:6{d) = 9Actci{d')- Note that 
LRAetei{d'; d) >0 and LRActciid; d) = 0, which yields oj{d; d) = 1. If 9Aete\{d') 
is close to 6Aete\{d)-, then LRAetei{d']d) is close to zero and uj{d';d) will be 
close to 1. However, if the difference between 9Actc\{d') and 6Aete\{d) is large, 
then LRActci{d';d) is large and uj{d']d) will be small. Thus, uj{d';d) defined 
in (2.15) truly characterizes the similarity between voxels d and d' . 

One substitutes g{zi{d),9{d);d) in (2.14) into (2.10) to build a new func- 
tion, denoted by l'Tcte\{9{d)]d), and then solves 0Tetei(c^) according to (2.11) 
by minimizing li(.^c\{9{d);d). Finally, to test Ho{d):R9{d) = bo, one uses 
g{zi{d),9{d);d) in (2.14) to calculate the TETEL ratio statistic LRT^ctciid) 
according to (2.13). Note that the key difference between LRTctci{d) and 
LR Actci{d) lies in their different sets of estimating equations. 

Although the two-stage procedure only combines the data in the voxels of 
N{d) with the data in voxel d, they may preserve the long-range correlation 
structure in the imaging data, because the neighborhoods of all voxels are 
consecutively connected. Thus, the two-stage procedure captures a substan- 
tial amount of spatial information in the imaging data. For the sake of space, 
we only present the asymptotic properties of 9Totci{d) and LRActciid) below. 

Theorem 2.2. If assumptions (A1)-(A3) and (A5)-(A7) in the supple- 
mentary document are true, then we have the following: 

(a) -v/ra(^Tetcl('^) — 9o{d)) converges to i'{d) = A^(0, S((i)) in distribution, 
where 9Q{d) is the true value of9{d) in the voxel d and Ti{d) = [D{d)V{d)~^ x 
D{d)'^]~^ , in which 

n 

D{d) = lim n-^y^de~g{zi{d)Md);d) 

i=l 
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and 

n 

V(d) = lim n-^V5(zi((i),eo(d);d)®'; 

n— j-oo ^— ' 
i=\ 

(b) under the null hypothesis HQ{d), LRictciid) converges in distribution 
to a x^(co) distribution. 

Theorem 2.2 establishes the asymptotic consistency and normality of 
^Tetei(c^) and the asymptotic distribution of Xi?Tetei(d)- Theorem 2.2 also 
shows that the asymptotic variance of 0Tetei(d) depends on all the data in 
N{d) U {d} for all subjects. Since the weights u}{d'; d) automatically put large 
weights on the neighboring voxels with similar pattern and small weights on 
the neighboring voxels with dissimilar pattern, it follows that the TETEL 
procedure produces more accurate parameter estimates and more powerful 
test statistics. 

TETEL has three features. TETEL not only downweights the data from 
the neighboring voxels with dissimilar signal pattern, but also incorporates 
the data from the neighboring voxels with similar signal pattern to adap- 
tively calculate parameter estimates and test statistics. TETEL avoids us- 
ing the same amount of smoothness throughout the whole image in most 
smoothing methods. Our theoretical results ensure the asymptotic consis- 
tency and normality of ^Tetci(d) and the asymptotic distribution of 
LRTctciid). Then, we can approximate the value of Z/iZ Tetei('^) at each 
voxel d. Finally, we correct for multiple comparisons by using either the 
family-wise error rate or false discovery rate (FDR) across the entire brain 
region [Benjamini and Hochberg (1995), Worsley et al. (2004)]. Since the 
smoothing stage in TETEL usually introduces the positive dependency 
among all LR-xctclid), it allows us to apply FDR in Benjamini and Yeku- 
tieli (2001) to control the false discovery rate. 

3. Simulation studies. Three sets of simulation studies were conducted to 
examine the finite sample performance of our AETEL and TETEL methods. 

3.1. Study I: Longitudinal data. We considered the following model: 

(3.1) Yij =/3o+ Pitij + (32Xi + PstijXi + bi + Sij 

for i = 1, . . . , n, where tij denotes time taking values in (1, 2, 3, 4, 5), Xi was 
independently generated from a A^(0, 1), and bi was independently generated 
from a A^(0,1). Errors Eij were independently generated from A^(0, 1) and 
X^(3) — 3, respectively, where X^(3) represents a chi-squared random variable 
with three degrees of freedom. The x^(3) — 3 distribution is very skewed and 
differs substantially from any symmetric distribution, such as a Gaussian 
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distribution. The true value of (/?o, /32)"^ was set at (1, 1, 1)"'" and (3^ was 
varied as 0, 0.05, 0.10, 0.15, and 0.20. We tested the hypothesis Hq: /33 = vs. 

: /33 / using LR Aetei- To assess both Type I and II error rates of LRActci, 
we used generahzed estimating equations assuming an exchangeable working 
correlation matrix to construct i-RActci and then compared it with the ETEL 
likelihood ratio statistic, denoted by Xi?Etch and the Wald statistic, denoted 
by Wn, obtained from the "true" linear mixed model (3.1) representing an 
ideal scenario. We considered n = 40, 60, and 80. The 1,000 replications 
were used to calculate the estimates of rejection rates with significance level 
a = 5%. 

The type I error rates of Xi?Actci and Wn are reasonably accurate for all 
sample sizes (n = 40, 60, or 80) considered and for all different distributions 
of error terms at the 5% significant level (Table 1). In contrast, the type I 
error rates of LR-^tei are slightly inflated for n = 40. The type II error rates 
for LRj^giei and Wn are similar under both error distributions and for all 
sample sizes (Table 1). However, the power of the three test statistics to 
reject the null hypothesis increases modestly when the distribution of the 
error terms follows the skewed distribution x^(3) — 3 (Table 1). This decline 
in the type II error rate was caused by the fact that the variance of x^(3) — 3 
is larger than that of A^(0, 1). Compared with Li?Actci and Wn, LR^^^i has 
slightly larger power, which may be due to its inflated type I error rates. 
Consistent with our expectation, the statistical power for rejecting the null 
hypothesis increases with the sample size n. 

3.2. Study II: Testing the type of time- dependent covariates. We used 
the simulation study for a type II time-dependent covariate in Section 4.1 
of Lai and Small (2007) to examine the finite sample performance of our 
AETEL method. The data were simulated under the mechanism 

yu = (3o + l^iXit + /32Xi,t-i + bi + eu and xu = fSsXi^t-i + ^it, 

where bi , eu , and eu are mutually independent and normally distributed 
with mean and variances 4, 1, and 1, respectively; the Xjf-process is sta- 
tionary, that is, XiQ ~ A^(0, (T^/(l — /3|)). We refer the reader to Lai and 
Small (2007) for more details. Note that xn is a type II covariate. We used 
our AETEL method with the following estimating equations: (a) the type 
II estimating equations according to (2.5), labeled type II; (b) the type III 
estimating equations according to (2.8), labeled type III; (c) GEE using the 
independent working correlation, labeled GEE independence; (d) GEE using 
the exchangeable working correlation, labeled GEE exchangeable; (e) GEE 
using the autoregressive AR-1 working correlation, labeled GEE AR-1. We 
compared the bias, root-mean-square error, and the efficiency of each case for 
the parameter /3i to the GEE independence case (the efficiency is the ratio 
of the mean-square error of the GEE independence case to that of the case). 
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Table 1 

Simulation study for comparing LRacicI, i-REtoi, an-d W„ for testing Ho: Pi =0 against 

Hi : A, / 











X^(3) -3 






iV(0,l) 




n = 40 


n = 60 


n = 80 


n = 40 


n = 60 


n = 80 




= 0.0 




0.078 


0.066 


0.059 


0.082 


0.070 


0.058 








0.066 


0.054 


0.055 


0.068 


0.064 


0.058 






w„ 


0.062 


0.064 


0.068 


0.078 


0.064 


0.054 




= 0.05 


LREtci 


0.112 


0.118 


0.118 


0.186 


0.280 


0.300 








0.088 


0.104 


0.102 


0.156 


0.254 


0.278 








0.094 


0.102 


0.100 


0.164 


0.244 


0.264 




= 0.10 


LREtci 


0.198 


0.182 


0.286 


0.548 


0.866 


0.804 








0.176 


0.160 


0.268 


0.506 


0.848 


0.792 






Wn 


0.168 


0.164 


0.270 


0.474 


0.786 


0.728 


/?3 


= 0.15 




0.280 


0.340 


0.394 


0.986 


0.930 


0.978 








0.250 


0.316 


0.374 


0.986 


0.920 


0.974 








0.268 


0.324 


0.356 


0.978 


0.892 


0.970 


/?3 


= 0.20 




0.560 


0.520 


0.720 


0.996 


1 


1 






LRActcl 


0.532 


0.488 


0.702 


0.990 


1 


1 






w„ 


0.512 


0.482 


0.682 


0.982 


0.998 


1 



Estimates of rejection rates were reported for A''(0, 1) and X'^(3) — 3 distributed data at 
3 different sample sizes (n = 40,60,80) at significance level a = 5%. For each case, 1,000 
simulated data sets were used. 



As we can see from Table 2, GEE independence and GEE AR-1 are biased, 
because they use some invalid estimating equations. The other three are all 
unbiased, with type II being more efficient than the other two. Combining all 
available valid estimating equations does improve efficiency. With the same 
type II estimating equations, our method has slightly less RMSE (0.0401 vs. 
0.0407) than Lai and Small (2007) 's method. 

Table 2 

Results of AETEL with various estimating equations for a type H time- dependent 

covariate 



Estimating equations Bias RMSE Efflciency 



Type II 


0.00 


0.040 


1.82 


Type III 


0.00 


0.053 


1.04 


GEE independence 


0.00 


0.054 


1.00 


GEE exchangeable 


-0.12 


0.104 




GEE AR-1 


-0.79 


0.661 
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Fig. 2. Two red regions of interest (ROIs) on a 30 x 30 image. The ROIs are indicated 
by the red area. 

3.3. Study III: Spatial data. We simulated data at all m = 900 pixels on 
a 30 X 30 phantom image (Figure 2). At a given voxel d, 

(3.2) yij{d) = /3o(d) + ^i{d)tij + ^2{d)xi + p^{d)tijXi + hi{d) + eij{d) 

for i = l,...,n and j = l,...,mj, where tij is the time taking values in 
(1,2,3,4,5), Xi was independently generated from a iV(0,l), and hi{d) was 
independently generated from a A^(0, 1). Errors £ij{d) were independently 
generated from A^(0, 1) and x^(3) — 3, respectively. We tested the hypotheses 
Hq : /33(d) = and Hi : /33(d) / across all pixels. To assess the Type I and 
II error rates at the pixel level, we set Po(d) = (3i(d) = /32(d) = across all 
pixels d and varied /33(d) as 0.0, 0.05, 0.10, 0.15, and 0.20. Specifically, we 
created two regions of interest (ROI) by setting /33(d) as 0.05, 0.10, 0.15, and 
0.20, and setting /33(d) = outside of the two ROIs in order to assess the 
finite sample performance of our method at different signal-to-noise ratios 
(SNRs). We considered n = 40 and 80. 

We used generalized estimation equations with an exchangeable working 
correlation matrix to calculate 6(d) and LRAetei(d) in Stage 1. In Stage 2 
we used the four first-order neighbors of pixel d to form N(d) and then 
calculated i-RTetei(rf)- As a comparison with the conventional analysis on 
image data, we first smoothed image data by using the heat kernel smoothing 
method with 16 iterations, which gave an effective smoothness of about 
4 pixels [Chung, Dalton and Davidson (2007)], and then calculated the Wald 
statistic based on GEE with an exchangeable working correlation matrix at 
each pixel. The 100 replications were used to approximate rejection rate 
with significance level a = 5%. 

As shown in Table 3, the Type I rejection rates outside of ROIs for both 
Li?Aetoi and Xi?xctci cire relatively accurate for all cases, while the statistical 
power for rejecting the null hypothesis in ROIs significantly increases with 
the absolute value of /33(d). Compared with Li^Aetoh -^-^Tctci has higher 
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Table 3 

Comparison of the two stages of TETEL for unsmoothed spatial data and the Wald test 
statistic for smoothed spatial data: true average rejection rates for voxels inside the ROI 
and false average rejection rates for voxels outside of the ROI were reported for N{0, 1) 
and X^(3) — 3 distributed data, and 2 different sample sizes (n = AO and 80) at a — 5%. 
For each case, 100 simulated data sets were used 



£i?Tetal Wald 











n 


40 


n = 


80 


n = 


40 


n = 


80 




Stags 




1 III. 


False 


True 


False 


True 


False 


True 


False 
















iV(0,l) 








0.05 


Staj 


?e 


1 


0.223 


0.088 


0.329 


0.068 


0.711 


0.101 


0.891 


0.105 




Staf 


le 


2 


0.302 


0.089 


0.426 


0.069 










0.10 


Staf 


le 


1 


0.571 


0.087 


0.820 


0.069 


0.964 


0.15 


0.991 


0.158 




Staf 


le 


2 


0.690 


0.088 


0.910 


0.070 










0.15 


Staf 


le 


1 


0.863 


0.089 


0.984 


0.069 


0.996 


0.184 


0.998 


0.177 




Staj 


re 


2 


0.954 


0.090 


0.998 


0.069 










0.20 


Staj 


le 


1 


0.987 


0.089 


0.999 


0.069 


0.999 


0.193 


0.999 


0.192 




Staf 


le 


2 


0.992 


0.090 


1.000 


0.069 


-3 








0.05 


Staf 


ye 


1 


0.117 


0.085 


0.122 


0.070 


0.313 


0.089 


0.331 


0.073 




Staf 


le 


2 


0.212 


0.090 


0.232 


0.070 










0.10 


Staf 


le 


1 


0.193 


0.087 


0.259 


0.069 


0.567 


0.099 


0.858 


0.099 




Staj 


le 


2 


0.278 


0.089 


0.411 


0.070 










0.15 


Staf 


re 


1 


0.313 


0.090 


0.447 


0.068 


0.847 


0.113 


0.948 


0.123 




Staj 


re 


2 


0.486 


0.091 


0.649 


0.070 










0.20 


Staf 


r,e 


1 


0.463 


0.090 


0.660 


0.069 


0.947 


0.130 


0.979 


0.145 




StaE 


re 


2 


0.653 


0.090 


0.859 


0.069 











statistical power for rejecting the null hypothesis in ROIs with /33(d) 7^ 0. 
In contrast, compared with Li?Actci and Li^xetei based on the unsmoothed 
imaging data, although the Wald statistic for the smoothed imaging data 
has higher statistical power for rejecting the null hypothesis in ROIs, its 
Type I error rate is inflated and increases with the absolute value of (33(d). 
The decline in the type I and II error rates is caused by the fact that the 
variance of x^(3) — 3 is larger than that of A^(0, 1). We also tried different 
degrees of smoothness and ROIs with different sizes and found that the 
degree of smoothness and ROI size can have profound effect on the Type I 
and II error rates of the Wald statistic (not presented here). 

4. Hippocampus shape. 

4.1. Hippocampus SPHARM-PDM representation. Let yjj((i) be the 3 x 
1 coordinate vector at voxel d on the left and right hippocampus SPHARM- 
PDMs and Xjj = (1, gender^, age^, SClj, SC2j, race Ij, race2j, iimejj)'^, where 
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SCI and SC2 were, respectively, dummy variables for haloper idol-treated SC 
patient and olanzapine-treated SC patient versus healthy controls, and racel 
and roce2 were, respectively, dummy variables for Caucasian and African 
American versus other race. Let yi{d) = {yaid)'^ , . . . ,yim^{d)'^)'^ and A(^B 
denote the Kronecker product of matrices A and B. We assume that the 
mean and covariance matrix of yi{d) are, respectively, given by 

E{y^{d)) = • • • m and Cov(yi,(d)) = Vi{d) = Ri{a{d)) ® 

where j5{d) is a 24 x 1 vector, Ri{a{d)) = (a(d)l-'~'^l) is the standard au- 
toregressive AR-1 correlation matrix and is a 3 x 3 covariance matrix 
of yij{d). We estimated a{d) and by using Pearson residuals, which 

were calculated by solving GEEs with an independent working correlation 
matrix. For now on, Vi{d) [or a{d) and are assumed to be known. For 

the data analysis, we used the moment model based on GEE in (2.2) since 
there is no time-dependent covariate except time itself. The g{zi{d),9{d);d) 
which is used in TETEL is given by 

n / xfi ® I3 \ ^ r / I3 

giz,{d),eidy,d) = Y,\ ■■■ y^idr' y^id)-l ■■■ 

Existing statistical methods of image data in SPM require that the error 
distribution is Gaussian and the variance is constant. The Shapiro- Wilk test 
rejects the normality assumption at many voxels of both the left and right 
hippocampus structures, and, thus, our nonparametric TETEL method is 
preferred for the analysis of this data set. 

Since our goal is to detect the difference in the SPHARM-PDM surface 
shape between the schizophrenia and control groups, we used LR^ctci and 
Li^Tctci to carry out the test. Moreover, in Stage 2, we used the closest neigh- 
bors of each voxel d to form N{d). The color-coded p- values of the Xi?Actci 
and Zi?xctci and their corrected p- values using FDR across the voxels of both 
the left and right reference hippocampi are shown in Figure 3 [Benjamini 
and Yekutieli (2001)], in which the top row is for the first stage (Xi?Actci) 
and the bottom row is for the second stage (ZiZxetei)- 

The analyses show strong shape differences in the superior, anterior parts 
of the left hippocampus, at the intersection of cornu ammonis 1 and cornu 
ammonis 2, previously not shown. Posterior shape changes at the hippocam- 
pal tail shown in chronic schizophrenics [Styner et al. (2004)] are detected 
here already in first episode patients. Furthermore, the results also confirm 
those reported in Narr et al. (2004) by indicating a strong medial shape 
difference in the central, left hippocampal body in first episode patients. 
Comparing the first and second rows, it is clear that TETEL shows advan- 
tages in detecting more significant and smoother activation 
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Fig. 3. Results from the longitudinal schizophrenia study. The first and third rows are 
for the first stage (LRacici): the color-coded raw p-value maps of group effect for the left 
hippocampus (a, b) and the right hippocampus (c, d), and the corresponding color-coded 
corrected p-value maps of group effect for the left hippocampus (i, j) and the right hippocam- 
pus (k, 1). The second and fourth rows are for the second stage (LRtcici): the color-coded 
p-value maps of group effect for the left hippocampus (e, f) and the right hippocampus (g, 
h), and the corresponding color-coded corrected p-value maps of group effect for the left 
hippocampus (m, n) and the right hippocampus (o, p). 



4.2. Hippocampus m-rep thickness. First, we considered the baseline anal- 
ysis. We used the moment model based on the estimating equations Xji(yji — 
where yn is the m-rep thickness measured at baseline for the ith sub- 
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(c) 




~i 1 1 1 r 
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Fig. 4. An m-rep model of a hippocampus: (a) an m-rep model of the hippocampus; 
(b) the boundary surface of the m-rep model of hippocampus; (d) m-rep radius ( or thick- 
ness) measures at the five atoms from two m-rep objects; (c) shows the — \ogjQ{p) -values 
for the Shapiro-Wilk test for the residuals at each atom on the left hippocampus; (e) shows 
the —logiQ{p) -values for the Shapiro-Wilk test for the residuals at each atom on the right 
hippocampus. The red horizontal line is the 0.05 cutoff line. 

ject at each medial atom of the left and right hippocampi; Xji is an 8 x 1 vec- 
tor given by Xji = (1, gender^, age^, SClj, SC2j, race Ij, race2 i,WBViif and 
/3 = (/3o,/3i, . . . j/Sy)'^. Existing statistical methods of image data in SPM re- 
quire that the error distribution is Gaussian and the variance is constant. 
The Shapiro-Wilk normality test was applied to check this parametric as- 
sumption of the general linear model at each atom for the left hippocampus 
and right hippocampus using the residuals. Figure 4(c) and (e) show that 
the Shapiro-Wilk test rejects the normality assumption at many atoms of 
both the left and right hippocampus structures, therefore, our nonparamet- 
ric AETEL method is preferred for the analysis of this data set. 

Here our goal is to detect differences in thickness of the hippocampus 
across the three groups. Hence, we set the null hypotheses Hq : /Ss = /34 = 
at all 24 atoms for both the left and right hippocampi. Accordingly, we have 

/O 1 0\ 

Vooooioooy 



EMPIRICAL LIKELIHOOD FOR LONGITUDINAL NEUROIMAGING DATA 21 
(a) Left (c) Left (e) Left 




5 10 15 20 5 10 15 20 0.0 0.1 0.2 0.3 0.4 

24 atoms 24 atoms pvmQIeft 



(b) Right (d) Right (f) Right 




1 1 1 — ^ 1 n 1 1 — ^ 1 I 1 1 1 1 

5 10 15 20 5 10 15 20 0.0 0.2 0.4 0.6 0.8 

24 atoms 24 atoms pvmQright 



Fig. 5. An m-rep model of a hippocampus: Maps of —logj^Q{p) -values for testing WBV 
as a type I time- dependent covariate (black) and a type II time-dependent covariate (red): 
(a) uncorrected —\ogj^Q{p) -values for left hippocampus; (b) uncorrected —\ogiQ{p) -values 
for right hippocampus; (c) corrected —logiQ{p) -values for left hippocampus; (d) cor- 
rected — \og^(j{p) -values for right hippocampus; (e) the goodness-of-fit test for the equation 
E{dp^i2{l3)\yiz — A*i3(/3)]} = for the 3rd atom on the left hippocampus; (f) the good- 
ness-of-fit test for the equation E{di3fii2{l3)['yi3 — l-iisiP)]} ~ /°'" 14f/i atom on the 
right hippocampus. 

and bo = (0,0)-^. We used LRxctci to carry out the test. The color-coded 
p-values of the LR^etci across the atoms of both the left and right refer- 
ence hippocampi are shown in Figure 5(a) and (b). The false discovery rate 
approach was used to correct for multiple comparisons, and the resulting 
adjusted p-values were shown in Figure 5(c) and (d). Before correcting for 
multiple comparisons, there was a significant group difference in m-rep thick- 
ness at the upper central atoms in the left hippocampus and some area in 
the right hippocampus. However, there is no significant group effect at any 
atom after correcting for multiple comparisons. 

Second, we did a longitudinal data analysis. The advantage of a longi- 
tudinal study over a baseline study is that it allows us to determine (i) 
whether the change patterns of the response are similar or not across the 
three groups; (ii) whether, on average over time, there is a difference in 
the response across the three groups. We considered the moment model 
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with Xjj = (1, gendeij, agej, SClj, SC2j, race Ij, race2j, WBVij,timeij, SClj * 
timeij,SC2i * timeij)'^ . 

Since the WBV is a time-dependent covariate, we need to verify its appro- 
priate type. Moreover, from a neuroscience point of view, the m-rep thickness 
at each atom serves as a local volumetric measure and covaries with WBV. 
We started with type III and used GEE in (2.2) with Vi = !{. Then we used 
the type II equations specified in (2.5) and tested whether WBV is type 
II against type III. The LRqy did not reject for almost all 24 atoms, sug- 
gesting WBV is a type II covariate for most atoms. Furthermore, we used 
the type I equations specified in (2.3) and tested whether WBV is type I 
against type II. The LRqp rejected that WBV was of type I for most atoms 
(Figure 5). This indicates the invalidity of some type I equations. We used 
goodness-of-fit statistics in Zhu et al. (2008b) to test whether some of the 
extra equations added for type I, such as 

E{di3if^isW)[yij -/iij(/3)]} = for aU s < j,j = 1, . . . ,mi, 

were not valid. For instance, for the 3rd atom on the left hippocampus, the p- 
value of the goodness-of-fit test for the newly added equation 
E{dj3^ijLi2{P)[yi3 — Att3(/3)]} = is smaller than 0.001 [Figure 5(e)]; for the 
14th atom on the right hippocampus, the p-value of the goodness-of-fit test 
for the newly added equation -E{i9/3,/ii2(/3)[yi3 — fJ-isif^)]} = is smaller than 
0.001 [Figure 5(f)]. Therefore, we treated WBV as a type II time-dependent 
covariate and used the corresponding estimating equations for the longitu- 
dinal data analysis. 

To determine whether the changing patterns of the thickness of the hip- 
pocampus over time are similar or not across the three groups, we tested 
the null hypotheses Hq : /Sg = /3iq = {(3g and /3io are the coefficients of the 
interaction terms of group and time) at all 24 atoms for each of the left hip- 
pocampus and the right hippocampus, and it turned out that the interaction 
terms were not significant for most atoms. Next we deleted the interaction 
terms and tried to look at whether there are differences in the responses 
across the three groups on average over time with respect to the null hy- 
potheses Hq : /Ss = /34 = at all 24 atoms for each of the left hippocampus 
and the right hippocampus. Again we only found that there was a signifi- 
cant difference through time in m-rep thickness at the upper central atoms 
in the left hippocampus across schizophrenia patients and healthy controls 
groups after correcting for multiple comparisons, but the differences were 
not significant at other atoms, nor at any atoms on the right hippocam- 
pus. The color-coded p- values of the LR A^tei across the atoms of both the 
left and right reference hippocampi are shown in Figure 5(e) and (f), and 
the corrected p- values are shown in Figure 6(g) and (h). Before correcting 
for multiple comparisons, there was a significant group difference in m-rep 
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Fig. 6. Results from the longitudinal schizophrenia study. The top row is for the baseline 
analysis: the color-coded uncorrected p-value maps of group effect for (a) the left hippocam- 
pus and (b) the right hippocampus; the color-coded corrected p-value maps of group effect 
for (c) the left hippocampus and (d) the right hippocampus after correcting for multiple 
comparisons. The bottom row is for the longitudinal analysis: the color-coded uncorrected 
p-value maps of group effect for (e) the left hippocampus and (f ) the right hippocampus; 
the color-coded corrected p-value maps of group effect for (g) the left hippocampus and (h) 
the right hippocampus after correcting for multiple comparisons. 

thickness at the upper central atoms in the left hippocampus, and the sig- 
nificance level is larger than that of the baseline analysis. Since the positive 
correlation is commonly observed in imaging data, we applied the false dis- 
covery rate (FDR) procedure in Benjamini and Yekutieli (2001) to correct 
for multiple comparisons. There is still a significant group effect at the upper 
central atoms in the left hippocampus. 

We compared the results by making the assumption that WBV was a type 
II time-dependent and also a type III time-dependent covariate. Treating 
WBV as a type II time-dependent covariate lowered the p-values, making 
some nonsignificant values for the group effect significant. On the other 
hand, we found that all the standard deviations associated with the param- 
eter estimates treating WBV as a type II time-dependent covariate were 
uniformly less than those treating WBV as a type III, which confirms that 
treating WBV as a type II gains efficiency by making use of more cor- 
rect estimating equations. Table 4 compares the standard deviations of the 
parameter estimates between treating WBV as a type II time-dependent 
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Table 4 

Standard deviation comparison of the parameter estimates between treating WBV as 
a type II time- dependent covariate and a type III time-dependent covariate at atom 11 of 

the left hippocampus 





Intercept 


Gender 


Age 


SCI 


SC2 


Racel 


Race2 


WBV 


Time 


Type III 
Type II 


0.367 
0.344 


0.078 
0.075 


0.007 
0.005 


0.062 
0.058 


0.058 
0.054 


0.097 
0.094 


0.102 
0.100 


0.237 
0.221 


0.022 
0.018 



covariate and a type III time-dependent covariate at atom 11 of the left 
hippocampus. 

The longitudinal analysis increased the significance level at those signifi- 
cant atoms for the group effect, compared to the baseline analysis. We were 
also able to observe the change difference across groups through time, al- 
though it is not much. Both the baseline analysis and longitudinal analysis 
suggest that there is an asymmetric aspect in that the left hippocampus 
shows larger regions of significance than the right one, and the significant 
positions of the group differences are around the lateral dentate gyrus and 
medial CA4 body regions for the left hippocampus. 

5. Discussion. We have developed TETEL for spatial analysis of neu- 
roimaging data from longitudinal studies. We have shown that AETEL al- 
lows us to efficiently analyze longitudinal data with different time-dependent 
covariate types. We have specifically combined all the data in the closest 
neighborhood of each voxel (or pixel) on a 3D volume (or 2D surface) with 
appropriate weights to calculate adaptive parameter estimates and adaptive 
test statistics. We have used simulation studies to examine the finite sam- 
ple performance of AETEL and TETEL. In our longitudinal schizophrenia 
study, we have used the boundary and medial shape of the hippocampus 
to detect differences in morphological changes of the hippocampus across 
time between schizophrenic patients and healthy subjects. For the m-rep 
thickness, we have found that WBV is an important time-dependent co- 
variate. Potential applications of our methodology include understanding 
normal and abnormal brain development, and identifying the neural bases 
of the pathophysiology and etiology of neurodegenerative and neuropsychi- 
atric disorders. 

Many issues still merit further research. One major issue is to develop 
a test procedure based perhaps on random field theory or resampling meth- 
ods to correct for multiple comparisons in order to control the family-wise 
error rate under the moment model (2.1). Another major issue is to extend 
the test procedure to conduct cluster size inference and examine its perfor- 
mance in controlling the Type I error rate. The test procedure may lead to 
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a simple cluster size test (cluster size test assesses significance for all sizes 
of the connected regions greater than a given primary threshold). Models 
with nonparametric components using TETEL also may prove to be useful 
directions to consider. 

SUPPLEMENTARY MATERIAL 

Proofs of Theorems 2.1 and 2.2 (DOL 10. 1214/11- AOAS480SUPP; .pdf). 
We present assumptions and proofs of Theorems 2.1 and 2.2. 

Acknowledgments. We thank the Editor, an Associate Editor, and two 
referees for valuable suggestions, which helped to improve our presentation 
greatly. We are thankful to Sylvain Gouttard, Steve Pizer, and Josh Levy 
for sharing their software. 

REFERENCES 

Almli, C. R., Rivkin, M. J., McKinstry, R. C. and Brain Development Cooper- 
ative Group (2007). The NIH MRI study of normal brain development (objective-2) ; 
Newborns, infants, toddlers, and preschoolers. Neurolmage 35 308-325. 

Arndt, S., Cohen, G., Alliger, R. J., Swayze, V. W. and Andreasen, N. C. (1991). 
Problems with ratio and proportion measures of imaged cerebral structures. Psychiatry 
Res. 40 79-89. 

Ashburner, J. and Friston, K. J. (2000). Voxel-based morphometry: The methods. 
Neurolmage 11 805-821. 

Beckmann, C. F., Jenkinson, M. and Smith, S. M. (2003). General multilevel linear 
modeling for group analysis in fMRI. Neurolmage 20 1052-1063. 

Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: A prac- 
tical and powerful approach to multiple testing. J. R. Stat. Sac. Ser. B 57 289-300. 
MR1325392 

BEN.IAMINI, Y. and Yekutieli, D. (2001). The control of the false discovery rate in 
multiple testing under dependency. Ann. Statist. 29 1165-1188. MR1869245 

Bowman, F. D., Caffo, B., Bassett, S. S. and Kilts, C. (2008). A Bayesian hierar- 
chical framework for spatial modeling of fMRI data. Neurolmage 39 146-156. 

Cao, J. and WORSLEY, K. J. (2001). Applications of random fields in human brain 
mapping. In Spatial Statistics: Methodological Aspects and Applications (M. MoORE, 
ed.). Lecture Notes in Statistics 159 170-182. Springer, New York. 

Chen, J., Variyath, A. M. and Abraham, B. (2008). Adjusted empirical hkelihood and 
its properties. J. Comput. Graph. Statist. 17 426-443. MR2439967 

Chung, M. K., Dalton, K. M. and Davidson, R. J. (2007). Tensor-based cortical surface 
morphometry via weighted spherical harmonic representation. IEEE Transactions on 
Medical Imaging 26 566-581. 

DiGGLE, P. J., Heagerty, P. J., LlANG, K.-Y. and Zeger, S. L. (2002). Analysis of 
Longitudinal Data, 2nd ed. Oxford Statistical Science Series 25. Oxford Univ. Press, 
Oxford. MR2049007 

Dryden, I. L. and Mardia, K. V. (1998). Statistical Shape Analysis. Wiley, Chichester. 
MR1646114 

DuvERNOY, H. (2005). The Human Hippocampus. Springer, New York. 



26 



X. SHI ET AL. 



Friston, K. J. (2007). Statistical Parametric Mapping: The Analysis of Functional Brain 

Images. Academic Press, London. 
Friston, K. J., Holmes, A. P., Poline, J. B., Price, C. J. and Frith, C. D. (1996). 

Detecting activations in PET and fMRI: Levels ol inlerence and power. Neurolmage 4 

223-235. 

Friston, K. ,J., Stephan, K. E., Lund, T. E., Morcom, A. and Kiedel, S. (2005). 
Mixed-effects and fMRI studies. Neurolmage 24 244-252. 

Hansen, L. P. (1982). Large sample properties of generalized method of moments esti- 
mators. Econometrica 50 1029-1054. MR0666123 

Hayasaka, S., Phan, L. K., Liberzon, L, Worsley, K. J. and Nichols, T. E. (2004). 
Nonstationary cluster-size inference with random field and permutation methods. Neu- 
rolmage 22 676-687. 

Hecke, W. v., Sijbers, J., Backer, S. D., Poot, D., Parizel, P. M. and Leemans, A. 

(2009). On the construction of a ground truth framework for evaluating voxel-based 

diffusion tensor MRI analysis methods. Neurolmage 46 692-707. 
Huettel, S. a., Song, A. W. and McCarthy, G. (2004). Functional Magnetic Reso- 
nance Imaging. Sinauer, Sunderland, MA. 
Imbens, G. W., Spady, R. H. and Johnson, P. (1998). Information-theoretic approaches 

to inference in moment condition models. Econometrica 66 333-357. MR1612246 
.Jones, D. K., Symms, D. K., Cercignani, M. and Howard, R. J. (2005). The effect 

of fiher size on VBM analyses of DT-MRI data. Neurolmage 26 546-554. 
Lai, T. L. and Small, D. (2007). Marginal regression analysis of longitudinal data with 

time-dependent covariates: A generalized method-of-moments approach. J. R. Stat. Soc. 

Ser. B Stat. Methodol. 69 79-99. MR2301501 
Liang, K. Y. and Zeger, S. L. (1986). Longitudinal data analysis using generalized 

linear models. Biometrika 73 13-22. MR0836430 

LlEBERMAN, J. A., TOLLEFSON, G. D., CHARLES, C, ZiPURSKY, R., ShARMA, T., 

Kahn, R. S., Keefe, R. S. E., Green, A. I., Gur, R. E., McEvoy, J., Perkins, D., 

Hamer, R. M., Gu, H. and Tohen, M. (2005). Antipsychotic drug effects on brain 

morphology in first-episode psychosis. Archives of General Psychiatry 62 361-370. 
LiNDQUlST, M. A. and Wager, T. D. (2008). Spatial smoothing in fMRI using prolate 

spheroidal wave functions. Hum. Brain Mapp. 29 1276-1287. 
Liu, Y. and Chen, ,J. (2010). Adjusted empirical likelihood with high-order precision. 

Ann. Statist. 38 1341-1362. MR2662345 
Logan, B. R. and Rowe, D. B. (2004). An evalution of thresholding techniques in fMRI 

analysis. Neurolmage 22 95-108. 
Luo, W.-L. and Nichols, T. E. (2003). Diagnosis and exploration of massively univariate 

neuroimaging models. Neurolmage 19 1014-1032. 
Luo, H. and Puthusserypady, S. (2005). A sparse Bayesian method for determination 

of flexible design matrix for fMRI data analysis. IEEE Trans. Circuits Syst. I. Regul. 

Pap. 52 2699-2706. 

Narr, K. L., Thompson, P. M., Szeszko, P., Robinson, D., Jang, S., Woods, R. P., 

Kim, S., Hayashi, K. M., Asunction, D., Toga, A. W. and Bilder, R. M. (2004). 

Regional specificity of hippocampal volume reductions in first-episode schizophrenia. 

Neurolmage 21 1563-1575. 
Newey, W. K. and Smith, R. J. (2004). Higher order properties of GMM and generalized 

empirical likelihood estimators. Econometrica 72 219-255. MR2031017 
Owen, A. B. (2001). Empirical Likelihood. Chapman & Hall/CRC, New York. 
Penny, W., Flandin, G. and Trujillo-Barreto, N. (2007). Bayesian comparison of 

spatially regularised general linear models. Hum. Bram Mapp. 28 275-293. 



EMPIRICAL LIKELIHOOD FOR LONGITUDINAL NEUROIMAGING DATA 27 



Pepe, M. S. and Anderson, G. L. (1994). A cautionary note on inference for marginal 
regression models with longitudinal data and general correlated response data. Comm. 
Statist. Simul. Comput. 23 939-951. 

Pizer, S. M., Fletcher, P. T., Joshi, S., Thall, A., Chen, J. Z., Fridman, Y., 
Fritsch, D. S., Gash, A. G., Glotzer, J. M., Jiroutek, M. R., Lu, C, 
MuLLER, K. E., Tracton, G., Yushkevich, p. and Chaney, E. L. (2003). De- 
formable m-reps for 3D medical image segmentation. Int. J. Comput. Vis. 55 85-106. 

POLINE, J. and Mazoyer, B. (1994). Analysis of individual brain activation maps us- 
ing hierarchical description and multiscale detection. IEEE Transactions in Medical 
Imaging 4 702-710. 

QlN, J. and Lawless, J. (1994). Empirical likelihood and general estimating equations. 

Ann. Statist. 22 300-325. MR1272085 
Qu, A., Lindsay, B. G. and Ll, B. (2000). Improving generalised estimating equations 

using quadratic inference functions. Biometrika 87 823-836. MR1813977 
Rogers, B. P., Morgan, V. L., Newton, A. T. and Gore, J. C. (2007). Assessing 

functional connectivity in the human brain by fMRI. Magn. Reson. Imaging 25 1347- 

1357. 

ROWE, D. B. (2005). Parameter estimation in the magnitude-only and complex- valued 

fMRI data models. Neurolmage 25 1124-1132. 
Salmond, C. H., Ashburner, J., Vargha-Khadem, F., Connelly, A., Gadian, D. G. 

and Friston, K. J. (2002). Distributional assumptions in voxel-based morphometry. 

Neurolmage 17 1027-1030. 
Schennach, S. M. (2007). Point estimation with exponentially tilted empirical likelihood. 

Ann. Statist. 35 634-672. MR2336862 
ShAfie, K., Sigal, B., Siegmund, D. and Worsley, K. J. (2003). Rotation space ran- 
dom fields with an application to fMRI data. Ann. Statist. 31 1732-1771. MR2036389 
Shi, X., Ibrahim, J. G., Lieberman, J., Styner, M., Li, Y. and Zhu, H. (2011). 

Supplement to "Two-stage empirical likelihood for longitudinal neuroimaging data." 

DOI:10.1214/11-AOAS480SUPP. 
Snook, L., Plewes, C. and Beaulieu, C. (2007). Voxel based versus region of interest 

analysis in diffusion tensor imaging of neurodevelopment. Neurolmage 34 243-252. 
Styner, M. and Gerig, G. (2003). Automatic and robust computation of 3d medial 

models incorporating object variability. Int. J. Comput. Vis. 55 107-122. 
Styner, M., Lieberman, J. A., Pantazis, D. and Gerig, G. (2004). Boundary and 

medial shape analysis of the hippocampus in schizophrenia. Med. Image Anal. 8 197- 

203. 

Styner, M., Lieberman, J. A., McClure, R. K., Weinberger, D. R., Jones, D. W. 
and Gerig, G. (2005). Morphometric analysis of lateral ventricles in schizophrenia and 
healthy controls regarding genetic and disease-specific factors. Proc. Natl. Acad. Sci. 
USA 102 4872-4877. 

Thompson, P. M., Cannon, T. D. and Toga, A. W. (2002). Mapping genetic influences 
on human brain structure. Annals of Medicine 24 523-536. 

Thompson, P. M. and Toga, A. W. (2002). A framework for computational anatomy. 
Comput. Vis. Sci. 5 13-34. 

WOOLRICH, M. W., Behrens, T. E. J., Beckmann, C. F., Jenkinson, M. and 
Smith, S. M. (2004). Multilevel linear modelling for fMRI group analysis using Bayesian 
inference. Neurolmage 21 1732-1747. 

Worsley, K. J., Taylor, J. E., Tomaiuolo, F. and Lerch, J. (2004). Unified univari- 
ate and multivariate random field theory. Neurolmage 23 189-195. 



28 



X. SHI ET AL. 



YuE, Y., LOH, J. M. and Lindquist, M. A. (2010). Adaptive spatial smoothing of fMRI 
images. Stat. Interface 3 3-13. MR2609707 

Zhu, H., Li, Y., Tang, N., Bansal, R., Hao, X., Weissman, M. M. and Peter- 
son, B. S. (2008a). Statistical modelling of brain morphological measures within family 
pedigrees. Statist. Sinica 18 1569-1591. MR2469324 

Zhu, H., Ibrahim, J. G., Tang, N. and Zhang, H. (2008b). Diagnostic measures for em- 
pirical likelihood of general estimating equations. Biometrika 95 489-507. MR2521595 

Zhu, H. T., Zhou, H., Chen, J., Li, Y., Styner, M. and Lieberman, J. (2009). Adjusted 
exponentially tilted likelihood with applications to brain morphology. Biometrics 65 



Biomedical Research Imaging Center 
University of North Carolina at Chapel Hill 
Chapel Hill, North Carolina 27599-7420 
USA 

E-MAIL: hzhu@bios.unc.cdu 



919-927. 



X. Shi 

J. G. Ibrahim 
M. Styner 
H. Zhu 

Department of Biostatistics and 



J. Lieberman 

New York State Psychiatric Institute 

10-51 Riverside Drive 

New York, New York 10032 

USA 



Y. Li 

St. Jude Children's Research Hospital 
262 Danny Thomas Place 
Memphis, Tennessee 38105-3678 
USA 



